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ABSTRACT 

We present an updated timing solution for the young, energetic pulsar 
PSR B1509— 58 based on 21.3 years of radio timing data and 7.6 years of X- 
ray timing data. No glitches have occurred in this time span, in contrast to other 
well-studied young pulsars, which show frequent glitches. We report a measure- 
ment of the third frequency derivative of v = (—1.28 ± 0.21) x 10^'^^ s^^. This 
value is 1.65 standard deviations from, i.e. consistent with, that predicted by the 
simple constant magnetic dipole model of pulsar spin-down. We measured the 
braking index to be n = 2.839 ±0.003 and show that it varies by 1.5% over 21.3 yr 
due to contamination from timing noise. Results of a low-resolution power spec- 
tral analysis of the significant noise apparent in the data yield a spectral index 
of a = —4.6 ± 1.0 for the red noise component. 

Subject headings: pulsars: individual (PSR B1509— 58), supernovae: individual 
(G320.4-1.2) 



1. Introduction 

Radio pulsars are powered by their rotational kinetic energy: as they emit electromag- 
netic radiation, their rotation rates decrease. Thus, measuring the rotational evolution of 
pulsars is a probe of the physics underlying these fascinating objects. Young, energetic 
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pulsars, such as PSR B1509— 58 are especially interesting due to their large spin-down lu- 
minosities and rapid spin-down rates. High rates of spin-down allow for the measurement of 
higher-order frequency derivatives, providing a deeper probe into the physics of spin-down. 
This slow-down can be described by 

z> = -Xi/", (1) 

where v = 1/P is the pulse frequency, z> is its derivative and n is the braking index. For 
braking by magnetic-dipole radiation in a vacuum, K is related to the dipole magnetic 
moment of the pulsar and n = 3 (Manchester & Taylor 1977). Differentiating (1) shows that 
n is given by 

. 2 

The braking index is only measurable for the youngest pulsars and values not entirely dom- 
inated by noise processes have been measured for only 5 objects. In all cases, the measured 
values are less than the canonical value of n = 3 (Lync ct al. 1993; Kaspi ct al. 1994a; Lyne 
et al. 1996; Deeter et al. 1999; Camilo ct al. 2000). Implications of observed braking indices 
less than 3 arc still being discussed (see for example, Blandford & Romani 1988) and include 
magnetic field growth or alignment with the spin axis, effects due to higher-order magnetic 
moments and the interaction of the pulsar with its relativistic particle wind. 

Most pulsars do not have measurable braking indices due to a combination of contami- 
nation from low-frequency noise processes and the relatively slow spin-down of middle-aged 
pulsars. Because of the form of the spin-down power law given in equation 1, successive fre- 
quency derivatives go approximately as powers of the first frequency derivative. Thus most 
pulsars, which have values of z> of ~ 10~^^s~^, should exhibit v on the order of ~ 10~^°s~^, 
which would require hundreds of years of observation to measure. Many younger pulsars 
that could conceivably have measurable v exhibit significant timing noise which masks the 
deterministic value of v due to dipole spin-down, resulting in measured values orders of mag- 
nitude larger than predicted, often with the wrong sign. Furthermore, these values are not 
stable and consequently do not accurately predict pulse arrival times. 

An important check on the validity of the spin-down law is the measurement of de- 
terministic higher-order frequency derivatives. An expression for ii is given by taking an 
additional derivative of equation (1), 

n(2n - l)z>3 



i/2 



(3) 



We define the second braking index to be mo = n(2n — 1). If equation 1 accurately charac- 
terizes the spin-down of the pulsar, then mo = m, where 



m = (4) 
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A measurement of m provides important insight into the spin-down of young pulsars and 
m 7^ mo imphes that the characteristic age is an under- or over- estimate of the true age, 
as well as having implications for a changing magnetic field (Blandford 1994). Evidently, 
the measurement of higher-order frequency derivatives will only be possible in the youngest 
pulsars that have large spin-down rates not dominated by noise processes. The Crab pulsar 
has a measured value of i) = (—6.45 ± 0.02) x 10^^^ s^"^ (Lyne et al. 1993), in agreement 
with the prediction from the spin-down law, in spite of probable contamination from timing 
noise and frequent glitches. 

PSR B1509— 58 has never been observed to glitch, hence is an excellent candidate to 
have a measurable u, as well as a very accurate measurement of n. Kaspi et al. (1994a) 
measured a value for i>' = (—1.02 ± 0.25) x 10"^^ s~^, based on 11 years of radio timing data, 
implying m — 14.5 ± 3.6, in agreement with the predicted value uiq — 13.26 ± 0.03. The 
uncertainty in i>' and thus m was known to be larger than the formal uncertainty due to 
contamination from significant timing noise in the data. To account, albeit not rigorously, 
for this effect, the quoted uncertainty is three times the formal uncertainty. 

In this paper, we extend the timing analysis performed by Kaspi et al. (1994a) using an 
additional 10.3 yr of radio data as well as 7.6 yr of X-ray timing data to increase the precision 
with which the frequency evolution of PSR B1509— 58 can be measured. This is accomplished 
by performing a phase-coherent analysis of the the data as well as a partially phase-coherent 
analysis that is less sensitive to the timing noise superimposed on the deterministic spin- 
down. In addition, we present a low-resolution spectral analysis of the timing noise and 
discuss possible physical causes of the noise process. 

2. Observations 

2.1. Radio Timing data 

In our 21-yr data span, observations were made with a variety of different data acqui- 
sition instruments, at two different radio telescopes. Observations of PSR B1509— 58 were 
conducted using the Molonglo Observatory Synthesis Telescope (MOST) at an observing 
frequency of 843 MHz from 1982 June 24 through 1988 June 23. Our analysis includes a 
total of 177 MOST pulse arrival times. For full details of these observations, see Manch- 
ester et al. (1985). All data used after 1990 were obtained using the Parkes 64-m telescope. 
Observations beginning 1990 March 15 and continuing at roughly monthly intervals until 21 
January 1994 were described in Kaspi et al. (1994a). The same observing system, consisting 
of a 64-channel analog filterbank spanning 320 MHz of bandwidth at a central frequency near 



-4- 



1400 MHz, was used until 1997 January. Our analysis includes a total of 127 pulse arrival 
times obtained with the 64-channel filterbank system at Parkes. From 1997 May through 
2003 October, filterbank observations were made near the same central frequency but using 
a 96-channel filterbank spanning 288 MHz of bandwidth. Details of this observing system, 
the central beam of the Parkes Multibeam receiver, can be found in Manchester et al. (2001). 
A total of 43 pulse arrival times were obtained with this system. Interspersed with the filter- 
bank arrival times are data obtained also at Parkes, but using the Caltech correlator system 
(Navarro 1994; Sandhu et al. 1997). Typical Parkes integration times were ~3000s for all 
systems. In spite of the variety of observing systems used, the Parkes data generally are of 
uniform quality, with typical folded pulse signal-to-noise ratios of ~20. 

The Parkes data set also includes a handful of observations at a significantly higher 
radio frequency. These data, obtained primarily at four epochs, were used to determine 
the dispersion measure (DM). Two observations at 2480MHz on 1990 March 16 were used 
with observations at 1400 MHz on 1990 March 15 and 1990 March 17 to determine the DM 
for the phase-coherent timing analysis and the DM remained constant within uncertainties 
throughout the entire data set (see Section 3.4). Folded pulse profiles were cross-correlated 
with a template to yield topocentric arrival times for each observation. Topocentric arrival 
times were fitted to a timing model (see Section 3.1) using the TEMPO ^ software package. 



2.2. X-ray timing data 

The results presented in the X-ray analysis were obtained with public data from the 
Proportional Counter Array (PGA; Jahoda et al. 1996) on board the Rossi X-ray Timing 
Explorer (RXTE). The PGA consists of an array of five collimated xenon/methane multi- 
anode proportional counter units (PGUs) operating in the 2-60 keV range, with a total 
effective area of approximately 6500 cm^ and a field of view of ~ 1° FWHM. We used 7.6 yr of 
archival RXTE observations collected in the "GoodXenonwithPropane" mode, which records 
the arrival time (with l-jis resolution) and energy (256 channel-resolution) of every unrejected 
xenon event as well as all the propane layer events. We used all xenon layers of each PGU 
in the 2-32 keV range because of the relative hardness of the source. The observations were 
reduced using software developed at MIT for handling raw spacecraft telemetry packet data. 
Data from the different PGUs were merged and binned at 1/1024 ms resolution. The data 
were then reduced to baryccntric dynamical time (TDB) at the solar system baryccntcr using 
the known position from radio intcrferometry (see Table 1) and JPL DE200 solar system 
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ephemeris. Each time series was folded with 64 phase bins using the contemporaneous radio 
timing ephemeris. Resulting pulse profiles were cross-correlated in the Fourier domain with 
a high signal-to-noise-ratio template created by adding phase-aligned profiles from previous 
observations. We implemented a Fourier domain filter by using only the first 6 harmonics 
in the cross-correlation due to the relatively sinusoidal pulse profile. The cross-correlation 
produces an average time of arrival (TOA) for each observation. 

3. Results 

3.1. Phase coherent timing analysis 

The stable rotation of pulsars allows (in the absence of glitches) a phase-coherent analy- 
sis of the timing data, that is, each turn of the pulsar is accounted for. This is accomplished 
by fitting the TO As to a Taylor expansion of pulse phase, (f). Pulse phase at any time t can 
be expressed as 

(Pit) = (Pito) + Mt - to) + ^Mt - to)^ + l^t - to)^ + ^1^'oit - to)". . ., (5) 

where the subscript '0' denotes a parameter evaluated at the reference epoch to- Both the 
radio and X-ray TOAs were fit to the above polynomial using the software package TEMPO. 
Figure 1 shows timing residuals from a phase-coherent timing analysis with TEMPO with 
radio timing data shown as dots and X-ray data shown as crosses. The position of the pulsar 
was held fixed at the position determined by radio interferometry (Table 1; Gaensler et al. 
1999). A constant phase offset was fitted between the radio and X-ray data to allow for 
the differences between the radio and X-ray pulse profiles. The top panel shows residuals 
with u, V. V and v fitted out. The middle panel shows the fourth frequency derivative also 
fitted; the bottom panel shows the residuals with the fifth frequency derivative fitted. The 
large concentration of power in the fifth frequency derivative is not predicted by theory, but 
is presumably due to timing noise. The absence of any sudden discontinuties visible in the 
residuals (Figure 1) indicates that no glitches have occured in the 21.3 yr of observations. 

The usual method of determining timing parameters involves fitting many higher-order 
derivatives to 'pre- whiten' the residuals (e.g. Kaspi et al. 1994b). However, we found that the 
value of V changes with each higher-order derivative fit, without converging to a single value, 
and without entirely 'whitening' the residuals. This is due to two effects: contamination from 
timing noise and covariance in the fitted parameters. The value of v is ambiguous from this 
type of analysis. A similar effect is seen with i>, though to a much smaller level (i.e. within 
formal uncertainties) and is thus unimportant. 
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Nevertheless, it is clear from the timing residuals that there is a very significant signal 
from V and it is unlikely that it is entirely due to timing noise. Hints that a frequency 
derivative may be dominated by a noise process are if the value is several orders of magnitude 
larger than the value predicted from theory, if it is of the wrong sign, and if it does not 
correctly predict pulse arrival times. None of these are true for PSR B1509— 58, indicating 
that though there is certainly some contamination from timing noise, it is likely that the 
noise component does not dominate. Thus to determine i>, we used a partially phase- 
coherent timing analysis as described in the next section. Spin parameters determined by 
the phase-coherent timing analysis are given in Table 1. 



3.2. Pcirtially phase- coherent timing analysis 

To determine v unambiguously, we employed a partially phase-coherent timing analysis. 
This method can also be used when a fully coherent solution is impossible due to glitches 
or large gaps in the data. We spht the 21.3-yr of data into subsets of approximately 2yr 
and fitted each phase-coherently. We ensured that the residuals were white for each interval 
and that the covariances between the fitted parameters {f, v and were less than 0.9. The 
results are shown in Figure 2. The slope of the line, determined by a weighted least squares 
fit, is v and has a value of (— 1.28±0.21) x 10^^^ s~^. This partially coherent method heavily 
samples the data points at the ends of each interval, while weakly sampling the middle, so 
we also considered intervals that better sampled the mid-point as much as possible. We used 
these latter data points only to confirm the measurement of v ^ but have not included them 
in the quoted value because it was impossible to optimize these points in the same manner 
as the first set. The uncertainty was obtained from a bootstrap analysis, as we suspected 
that the formal uncertainty underestimated the true uncertainty due to contamination from 
timing noise. The bootstrap is a robust method of determining errors when a small number 
of sample points is available, as in this case (Efron 1979). 

In addition to measuring i>, we measured the braking index for each subset. The results 
are shown in Figure 3. The weighted average value is 2.839 ±0.003. A weighted hnear least- 
squares fit showed that ii is consistent with zero, i.e. that there is no underlying constant 
trend in the data. There is, however, significant variation from the average value; the reduced 
is 15 for 9 degrees of freedom. This value of n is consistent with both the measurement 
from the fully coherent timing analysis (see Table 1) and that measured by Kaspi et al. 
(1994a). 
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3.3. Timing noise 

Despite the renowned steady rotation of pulsars, there is timing noise apparent in most 
pulsar timing data. A low-frequency timing noise process can be seen in the residuals of 
PSR B1509— 58 (Figure 1). These processes are typified by long-term polynomial-like trends 
in timing residuals after all determinstic spin-down effects have been removed. Timing noise 
has been shown to be correlated with P, thus tends to be most apparent in young pulsars 
(Cordes & Downs 1985; Arzoumanian et al. 1994). The low- frequency nature of the noise 
process (giving rise to the adjective 'red'), complicates both its removal and understanding. 

One common way of characterizing timing noise is with the Ag parameter, given by 
Ag —log{-^\i>\t^) (Arzoumanian et al. 1994). This parameter estimates how much the timing 
noise in the i) term contributes to the cumulative phase of the pulsar, assuming that the 
measurement of i) is dominated by timing noise. One must be careful when calculating 
the Ag parameter for pulsars that have deterministic values of i), since these are no longer 
dominated by a noise process. Instead, the contribution to i) from timing noise must be found. 
We estimated the noise in this term using the partially coherent analysis explained in the 
previous section, by calculating the root mean square of the scatter about the deterministic 
trend. However, we found the uncertainty in the rms value to be larger than the rms value 
itself, so the value of Ag is consistent with zero. This is not inconsistent with the prediction 
made by Arzoumanian et al. (1994), given the large amount of scatter in their data. 

For a general review of timing noise in pulsars, see papers by Cordes (1980), Cordes 
& Helfand (1980) and Cordes & Greenstein (1981). The physical causes of timing noise in 
pulsars are poorly understood; possibilities range from changes in moment of inertia due to 
random pinning and unpinning of superfluid vortices in the core of the neutron star (Cordes 
& Greenstein 1981) to torques acting on the crust due to interactions between the pulsar and 
its magnetosphere (Cheng 1987), or that isolated pulsars may be experiencing free precession 
(Stairs et al. 2000). 

Another possible cause of apparent timing noise in some objects is the presence of one 
or more planet-mass objects orbiting the neutron star. If these planets are of sufficiently 
low mass, their presence may appear to be 'timing noise' for many years. In fact, this 
possibility has been suggested by Rots (2002, 2004) for PSR B1509— 58. He performed a 
fit consisting of five orbits of ^solar mass planets to the same RXTE timing residuals that 
we have used in our analysis. To determine if these orbits were stable, we performed a 
periodogram analysis on both the X-ray timing residuals and the entire data set. We found 
roughly similar periodicities in the 7.6-yr X-ray data set, but these periodicities are not seen 
in the entire 21.3-yr data set. This indicates that the 'periodicities' are not stable in time 
and thus hkely do not represent orbits of planets around PSR B1509— 58. 
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3.4. Dispersion Measure Variations 

One possible contributing factor to the observed timing noise is changes in the DM of 
the pulsar. Variations in DM have been measured for several pulsars, notably for the Crab 
and Vela pulsars (Backer et al. 1993; Hamilton et al. 1977). It has been postulated that the 
short-term DM variations observed in the Crab are due to wisps passing through the line 
of sight. A typical variation in DM for the Crab pulsar is ADM=0.02 pc cm~^ with a rise 
time of ~ 40 days and decay time of ~ 500 days, although more recently, a DM variation of 
ADM=0.15pc cm~^ was reported (Wong et al. 2001). This caused a change in arrival time 
of ~ 1ms, at an analysis frequency of 327 MHz. A similar change in DM in PSR B1509— 58 
would result in a change in arrival time of ~ 0.32 ms at 1400 MHz, undetectable in our data. 

The DM for PSR B1509— 58 was determined at four epochs where sufficient multifre- 
quency data were obtained as shown in Table 2. Each value was determined by fitting only 
for V and DM using only the multifrequency data occuring over a short time span (see Table 
2). The DM was held fixed for the remainder of the timing analysis at the value obtained for 
the 1990 epoch. Unfortunately, any changes in the DM that occured at epochs other than 
those with multifrequency data remain unmodelled and will thus contribute to the noise in 
the data, since the time of arrival of a pulse changes with DM. 

A weighted least-squares fit was performed on the four DM measurements and the best- 
fit slope gives ADM= (0.42 ± 0.19) pc cm~^yr~^. Data spanning these four epochs were 
also fitted phase-coherently using TEMPO resulting in ADM = (0.76 ± 0.25) pc cm~^yr~-'^. 
Although this second measurement appears marginally significant, the small variations in 
observing frequency throughout the data allow timing noise to be absorbed into the ADM 
measurement, thus underestimating the uncertainties. Prom Table 2, we can see that the 
largest single change in DM (occuring between 1990 and 1996), ~ 4pc cm~^, would cause a 
delay in arrival time of ~ 8.5 ms at 1400 MHz. There is no evidence, however, of a sudden 
change in DM, and a gradual change in DM of this magnitude would be difficult to detect 
due to the instrinsic scatter in the arrival times of the pulses. 

Another determination of DM is obtained for the time period starting in 1996, when 
RXTE began observing the source. X-rays are unaffected by changes in DM, thus by com- 
paring the X-ray and radio residuals, we can determine if changes in DM have occurred. 
We fit the X-ray data and all contemporaneous radio data with TEMPO, fitting out many 
higher-order derivatives to 'pre-whiten' the residuals. Though it is impossible to fit for DM 
directly with TEMPO between radio and X-ray data, we can compare the residuals. Due to the 
uncertainty between the radio and X-ray pulse profiles, we are required to fit a constant offset 
between the two data sets. Thus we are insensitive to any constant change in DM between 
this and an earlier epoch. The timing residuals show no evidence for a systematic change 
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in DM over 7.6 yr. However, the scatter remaining in the data is larger than expected given 
the TOA uncertainties. This scatter may be, in part, due to changes in DM on timescales 
of ~l-2 months, appearing random in nature in our sparsely sampled data set. The paucity 
of multifrequency data results in the inability to discern between short-term DM variations 
and a random noise process that is independent of observing frequency. We can rule out 
changes in DM that would change the arrival time of the pulses by more than the scatter 
in the data, ~10ms. Thus an observing frequency of 1400 MHz, DM variations larger than 
4.7 pc cm"^ are ruled out. 



Due to the low-frequency nature of the noise process active in most pulsars, as well as 
the uneven sampling of our data, a traditional Fourier spectral analysis is rendered ineffective 
as an estimator of the power spectrum (Deeter & Boynton 1982). These problems can be 
partially circumvented by using a method of extracting power spectra based on orthogonal 
polynomials estimators developed by Groth (1975) and Deeter (1984). This method provides 
a low-resolution but reasonable estimate of the power contained in timing noise over a range 
of frequencies. To extract a power spectrum from the noise component of the timing data, 
the deterministic spin-down terms (in this case, 0, z/, z>, i) and v) are fit out with TEMPO. The 
next order term contains a large contribution from red noise (Figure 1), while higher-order 
terms will contain declining contributions from red noise and a white noise contribution 
attributable to measurement errors. 

First, it is necessary to construct orthonormal polynomials pj,j — 0, 1, .., such that 



where each ti represents one TOA. Due to the nature of complete sets of orthogonal polyno- 
mials, they can, given the proper coefficients, describe any function. In this case, they are 
used to describe the timing residuals: 



3.5. Power Spectrum Analysis 




(6) 




(7) 



where R(t) are the residuals, and the coefficients are fitted for and described by: 




(8) 



Since 0, i/, z>, i> and i>, have been already been fitted, they are completely covariant with 
Co,Ci,C2,C3 and C4 and thus contribute no additional information about the noise process. 



-10- 



Therefore C5 is the lowest order coefficient tliat contains useful information about the timing 
noise. An estimate of the power contained at the lowest frequencies, that is / = 1/T, where 
T is the length of the data span, is given by Sm = jcsp (Groth 1975). It is possible to obtain 
further estimates of the power at higher frequencies by splitting the data into m = 2, 4, 8. . . 
sections, which give estimates of the power contained at / = m/T. For each value of m, 
there are m estimates of Sm, which are averaged to give the mean power estimate. Each 
average of S^. is distributed as a distribution with m degrees of freedom. The averages 
are then corrected by a factor to account for the difference between the median and the mean 
of a distribution (Kaspi et al. 1994b). The results of this analysis are plotted in Figure 
4. The highest frequency term is certainly due to measurement errors, and is of the same 
approximate power level as the previous three terms. Thus, these terms are not included 
in the fit of the power spectrum. The slope obtained by fitting the three lowest frequency 
terms is a = —4.6 ± 1.8. However, fitting the fourth term, (whose uncertainty places it 
just above the white noise level) gives the same spectral index while better constraining the 
uncertainties, and implies a spectral index oi a = —4.6 ± 1.0. 



4. Discussion 

4.1. Timing Analysis 

The measurement of the third frequency derivative, i>' = (—1.28 ± 0.21) x 10~^^s~^ is 
1.65 standard deviations from the value predicted by the spin-down law. Hence the observed 
value and the model prediction are consistent. Typically, uncertainties on i>' fall as 
though the prominence of timing noise in this case will undoubtably lengthen the amount 
of time required. Based on the bootstrap analysis used to determine uncertainties, another 
~ 10 yr of observations are required to reduce uncertainties by a factor of ~ 2. Assuming 
that we have measured a stable value of i) and that the uncertainties will be reduced by a 
factor of 2, then the simple spin-down law derived assuming a constant magnetic field may 
not sufficiently describe the spin-down of PSR B1509— 58. 

The consequences of a measured i) different from the predicted value are most clear when 
comparing the measured and predicted values of the second braking index. The measured 
value of ]) implies a second braking index m = 18.3 ±2.9, larger than tjiq = 13.26 ±0.03, but 
not significantly so given the uncertainties. Blandford (1994) discusses the implications of 
m 7^ mo, namely that X is a function of time. This discussion is particularly useful because 
it does not require knowledge of the true braking index, uq, which may be different from 
the measured value and cannot be predicted from theory. However, since we do not make 
any assumption about no, we cannot say whether K is changing at the present time. A 
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measurement of m 7^ mo does give us insight into the variation of K in the pulsar's past. In 
particular, a value of m < mo implies a period of magnetic field growth early in the pulsar's 
history and a true age larger than the characteristic age, while m > mo reduces the inferred 
age of the pulsar. Given m — mo = 5.04 ± 2.99, and using Blandford's prescription for the 
inferred age of the pulsar, the true age of PSR B1509— 58 is 1000±700yr, just over half of 
the characteristic age, Tc = 1691 yr. 

The supernova remnant (SNR) believed to be associated with PSR B1509— 58 was origi- 
nally estimated to have an age of 6 — 21 kyr (Seward ct al. 1983), indicating that the pulsar's 
characteristic age was underestimating the age of the system by a factor of ~ 3.5 — 12. 
However, Gaensler et al. (1999) showed that the SNR may be as young as the characteristic 
age of the pulsar. They show that the SNR and pulsar are interacting and explain the large 
size and unusual morphology of the remnant by hypothesizing that the remnant is due to an 
explosion of high kinetic energy or low ejected mass occurring on the edge of a cavity, thus 
accelerating the expansion of the remnant. At first glance, the measured value of m appears 
to worsen the possible age discrepancy between the pulsar and the SNR. However, the age 
of the remnant is clearly not well determined and could be as low as 1000 yr. 

Blandford & Romani (1988) discuss the imphcations of a constant observed braking 
index, n, where in this case we must make some a priori assumption of the true braking 
index, typically assumed to be ng = 3. They give several reasons why n < 3 may be 
observed. The torque may scale with frequency differently, either due to the outflow of 
plasma removing a signiflcant amount of angular momentum, or a non-dipolar magnetic 
field. Another interpretation of n < 3 is that the magnetic moment is currently varying 
in time, that is, K — K{t). If this is the case, we can then determine a time scale for 
torque evolution and a check on the functional form of evolution by finding the first two 
dimensionless derivatives ol v — —K{t)v'^: 

oi = n - no = — -, (9) 
A V 

and 



(i2 = m - no[2no - 1 + 3(n - no)] = — ^. (10) 

Equation (9) gives a time scale for torque evolution while equation (10) provides a check for an 
assumed functional form for K{t). For PSR B1509-58, di = -0.161±0.003 and = 3.5±3. 
For exponential growth of this measurement of d\ implies a time scale for growth of 
~ 20 000 yr. If the form of K is exponential we would expect that ^2 = (c^i)^ = 0.026, 
which cannot be ruled out with the present uncertainties for d-i- We find that the present 
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measurements are similarily in agreement with a power-law form for K{t). The measurement 
of m and the large uncertainties, however, are consistent with the simplest explanation of a 
constant value of K and thus a constant magnetic field. 

Despite scatter due to timing noise, the braking index is constant over 21.3 yr to within 
1.5%. This result is similar to that of the Crab pulsar's braking index, reported to be 
constant to within 0.5%, where measurements were obtained between glitch events and are 
typically based on five years of data (Lyne ct al. 1993). This reinforces the need for long-term 
timing to measure the true value of braking indices to such precision. 

4.2. Timing noise analysis 

4.2.1. Glitches 

Our phase-coherent analysis shows that, due to the absence of any sudden discontinuities 
in the timing residuals, PSR B1509— 58 has not glitched in 21.3 yr. PSR B1509— 58 is the only 
young pulsar known that has never glitched over a long period of time. Lyne (1999) showed 
that there is a correlation not only between P and timing noise, but P and glitch activity, 
although some of the youngest pulsars (i.e. the Crab) have glitch activity parameters smaller 
than expected from the correlation. The explanation for this is that the internal temperature 
of the youngest pulsars is too high for large scale vortex pinning and unpinning. If this is 
the case, then the fact that PSR B1509— 58 has the lowest glitch rate among young pulsars 
would suggest that it has the highest internal temperature, in spite its characteristic age 
being larger than that of the Crab pulsar. On the other hand, so-called 'anomalous X-ray 
pulsars' (AXPs) have been observed to glitch with parameters that are not very different from 
those seen in radio pulsars (Kaspi et al. 2000; Kaspi & Gavriil 2003; Dall'Osso et al. 2003). 
Given their estimated surface temperatures and luminosities (measured via X-ray spectral 
observations) AXPs appear to be much hotter than any known radio pulsar. If AXP glitches 
indeed have the same physical origin as do glitches in rotation-powered pulsars, this argues 
that the glitch phenomenon is not strongly affected by neutron-star temperature. 

4-2.2. Spectral Analysis 

The spectral index for the timing noise of PSR B1509— 58, a — — 4.6 ± 1.0 is marginally 
steeper than that measured for other pulsars. The low-resolution nature of the Deeter 
method, however, provides only a rough estimate of the true spectral index of the noise. The 
spectral index of the timing noise for the Crab pulsar has been determined several times. 
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the most recent of which found a = — 3 as well as a periodic component with period of 
568 ± 10 days (Scott et al. 2003), using a modified Fourier analysis (possible in this case due 
to the dense sampling of the Crab timing data from 1982-1989). The spectral index was 
then verified with the Deeter polynomial method, though this analysis does not probe the 
lowest frequencies (i.e. longest timescales) available. Previous analyses of Crab timing data 
were consistent with a ~ — 4 (Deeter 1981; Cordes 1980). Cordes & Helfand (1980) show 
that the timing noise of 11 pulsars is consistent with a pure random walk in phase, frequency 
or frequency derivative. However, D'Alessandro et al. (1995) show in a study of 45 pulsars 
that only 5 are consistent with a pure random walk, while most pulsars exhibit timing noise 
either due to a combination of a random walk and discrete jumps, or inconsistent with either 
hypothesis. Other pulsars have been shown to have spectral indices of timing noise ranging 
from OL ~ —0.5 to —2.4 (Baykal ct al. 1999). Thus there is a range of possible spectral indices 
that describe timing noise in different pulsars. 

Although there is no physical reason why the spectral index must be an even integer, 
these are the only cases for which analytical solutions exist. For PSR B1509— 58, the mea- 
sured spectral index is consistent with having a component of frequency noise, that is, a 
random walk in pulse frequency corresponding to a = —4. The measured value of the spec- 
tral index indicates that there may also be a component due to a random walk in frequency 
derivative, characterized by a = —6. Although the measured spectral index alone would 
seem to indicate a random walk in frequency, the fact that the residuals are not Gaussian 
distributed after fitting for higher derivatives indicates that this is not the case. A random 
walk in the /cth derivative of phase is mathematically equivalent to a white noise process in 
the (fc -|- 1) derivative (Groth 1975). This should correspond to 'white' residuals after fitting 
out the [k -|- 1) derivative, not seen in the timing residuals of PSR B1509— 58. 

One of the problems in quantifying timing noise is the lack of a comprehensive theory of 
its physical causes. In fact it is possible that there are several causes of timing noise active 
in one object, thus complicating the interpretation of the measured spectral index. Free 
precession is an interesting solution to the timing noise problem because it can be detected 
from the changes in the pulse profile as well as a periodicity in the timing residuals (Stairs 
et al. 2000). Unfortunately, the data quality and lack of a long baseline for many pulsars 
prohibits the discovery of this process. If undetected, precession would contribute to the 
apparent timing noise of the pulsar and the spectral index of the noise, despite the fact that 
it is not a random process. Additionally, other physical causes of noise, (e.g. magnetospheric 
torques or random vortex pinning and unpinning) will contribute to the measured spectral 
index. 

Another common problem in assessing the cause of timing noise is that quasi-periodicities 
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are apparent in the timing residuals of most pulsars. An artificial quasi-periodicity will 
change with the length of the data set, while a real quasi-periodicty will not (Konacki et al. 
1999). To begin to truly understand timing noise in pulsars, it is crucial that new high- 
resolution spectral methods, capable of discerning discrete features due to periodicities, are 
developed. New methods must be able to deal with both the spectral leakage and data sam- 
pling problems that presently constrain our measurements and thus understanding of timing 
noise. Another problem in developing high-resolution methods for very long period noise is 
the extremely long time baselines required to observe many periods of the noise; the only 
way to solve this is by brute force, that is, by timing pulsars for hundreds of years. A better 
understanding, both mathematically and physically, of timing noise will allow for a better 
understanding of the deterministic spin-down of pulsars underlying the noise processes. 
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Fig. 1. — Timing residuals for PSR B1509— 58. Radio TOAs are shown as dots; the X-ray 
TOAs as crosses. The top panel has a quartic polynomial (i.e. z/'), removed, the middle 
panel has a quintic removed and the bottom panel shows the residuals after the removal of 
a sixth degree polynomial. See the electronic edition of the Journal for a colour version of 
this figure. 
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Fig. 2. — Second derivative of pulse frequency from phase-coherent subsets versus epoch. 
The slope of the line is v and has a value of (-1.28 ± 0.21) x lO^^i s""^. 
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Fig. 3. — Braking index calculated at each epoch. There is no statistically significant change 
over 21.3 years of data. The average value is 2.839 ±0.003, in agreement with the previously 
reported value (Kaspi ct al. 1994a) and the value obtained in a fully phase-coherent analysis 
(Table 1). The reduced value is 15 for 9 degrees of freedom, suggesting contamination by 
timing noise (see Section 3.2). 



-20- 




Fig. 4. — Power spectrum of the timing noise exhibited by PSR B1509— 58. The best-fit hne 
gives a spectral index ofa = — 4.6±1.0. The last three points are dominated by measurement 
errors, thus are not fit as part of the red noise spectrum. 
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Table 1. Parameters for PSR B1509-58. 



Parameters for phase- coherent analysis. 



Dates (Modified Julian Day) 


45114 - 52925 


Epoch (Modified Julian Day) 


49034.5 


Right Ascension ^ (J2000) 


15h 13m 55.62s 


Declination ^ (J2000) 


-59°08'09.0" 


V (Hz) 


6.633598804(3) 


z> (10-11 s-2) 


-6.75801754(4) 


i) (10-21 s-3) 


1.95671(2) 


Braking Index, n 


2.84209(3) 


Dispersion Measure ^ (pc cm-^) 


253.2 


Parameters for partially coherent analysis. 


Braking Index, n 


2.839(3) 


v (10-31 s-4) 


-1.28(21) 


Second braking index, m 


18.3(2.9) 


'^Held fixed for phase-coherent analysis. 


Table 2. Dispersion Measure at four epochs. 


Epoch 


DM (pc cm-3) 


March 15 1990 


253.2±1.9 


February 1 - February 4 1996 


257.4±1.2 


February 29 - March 1 1996 


255.5±0.9 


January 12 - January 14 1997 


254.8±1.1 



